/* yyt program's support - additional functions, classes
 * * Gutzwiller factors
 * * Distributions finction
 * * Dispersion function
 */
#include"yyt05Funcs.h"
using namespace std;

/* Functions from Yuan, Yuan, Ting paper */

yytFunctions::yytFunctions(Parameters *par)
{
	setPar(par);
}

yytFunctions::yytFunctions(Parameters *par, Parameters *X0)
{
	setPar(par);
	setX0(X0);
}

int yytFunctions::setPar(Parameters *par)
{
	if(par->dim() < 5) { return -1; }
	sdelta = (*par)[0]; // 0. sdelta = 1-n (small delta),
	J = (*par)[1]; // 1. J (exchange coupling constant)
	t = (*par)[2]; // 2. t (hopping integral for the nearest neighbors),
	U = (*par)[3]; // 3. U (repulsive potential for electrons in the same spatial point).
	beta = (*par)[4]; // 4. beta (inverse temperature),
	L = (*par)[5]; // 5. L (siza of lattice)

	n = 1.-sdelta;
	return 0;
}

int yytFunctions::setX0(Parameters *X0)
{
	if(X0->dim() < 5) { return -1; }
	BDelta = (*X0)[0];
	chi = (*X0)[1];
	m = (*X0)[2];
	d = (*X0)[3];
	mu = (*X0)[4];

	r = n/2. + m;
	w = n/2. - m;
	return 0;
}

double yytFunctions::gt()
{
	double part1 = (n-2.*d)/(n-2.*r*w);
	double part2 = sqrt((1.-w)*(1.-n+d)/(1.-r)) + sqrt(w*d/r);
	double part3 = sqrt((1.-r)*(1.-n+d)/(1.-w)) + sqrt(r*d/w);
	//fpri2wntf(stderr, ":: %f, %f, %f, %f\n", r, w, n, d);
	return part1*part2*part3;
}

double yytFunctions::gs()
{
	return pow((n-2.*d)/(n-2.*w*r),2.);
	//double p1 = (n-2.*d)/(n-2.*r*w);
	//return p1*p1;
}

double yytFunctions::dgtdm()
{
	return -(((sqrt(((1 + d - n)*(1 + m - n/2.))/(1 - m - n/2.)) +
	          sqrt((d*(-m + n/2.))/(m + n/2.)))*(-2*(-m + n/2.) + 2*(m + n/2.))*
	        (sqrt(((1 + d - n)*(1 - m - n/2.))/(1 + m - n/2.)) +
	          sqrt((d*(m + n/2.))/(-m + n/2.)))*(-2*d + n))/
	      pow(-2*(-m + n/2.)*(m + n/2.) + n,2)) +
	   ((sqrt(((1 + d - n)*(1 + m - n/2.))/(1 - m - n/2.)) +
	        sqrt((d*(-m + n/2.))/(m + n/2.)))*
	      ((-(((1 + d - n)*(1 - m - n/2.))/pow(1 + m - n/2.,2)) -
	           (1 + d - n)/(1 + m - n/2.))/
	         (2.*sqrt(((1 + d - n)*(1 - m - n/2.))/(1 + m - n/2.))) +
	        (d/(-m + n/2.) + (d*(m + n/2.))/pow(-m + n/2.,2))/
	         (2.*sqrt((d*(m + n/2.))/(-m + n/2.))))*(-2*d + n))/
	    (-2*(-m + n/2.)*(m + n/2.) + n) +
	   ((((1 + d - n)/(1 - m - n/2.) +
	           ((1 + d - n)*(1 + m - n/2.))/pow(1 - m - n/2.,2))/
	         (2.*sqrt(((1 + d - n)*(1 + m - n/2.))/(1 - m - n/2.))) +
	        (-((d*(-m + n/2.))/pow(m + n/2.,2)) - d/(m + n/2.))/
	         (2.*sqrt((d*(-m + n/2.))/(m + n/2.))))*
	      (sqrt(((1 + d - n)*(1 - m - n/2.))/(1 + m - n/2.)) +
	        sqrt((d*(m + n/2.))/(-m + n/2.)))*(-2*d + n))/(-2*(-m + n/2.)*(m + n/2.) + n);

//	double sqrt1a = sqrt((1.-w)*(1.-n+d)/(1.-r));
//	double sqrt1b = sqrt(w*d/r);
//	double sqrt1 = sqrt1a + sqrt1b;
//	double sqrt2a = sqrt((1.-r)*(1.-n+d)/(1.-w));
//	double sqrt2b = sqrt(r*d/w);
//	double sqrt2 = sqrt2a + sqrt2b;

//	double c =  dgsdm()*gt()/(2.*gs()) +
//			0.5*gt()/sqrt2*(sqrt1a*(1./(1.-w)+1./(1.-r)) - sqrt1b*(1./w+1./r)) +
//			0.5*gt()/sqrt1*(-sqrt2a*(1./(1.-w)+1./(1.-r)) + sqrt2b*(1./w+1./r));
	//printf("............. %e, %e OK\n", aa-c, aa-b);

//	return 0.5*gt()*(dgsdm()/gs() + (1./(1.-w)+1./(1.-r))*(sqrt1a/sqrt2 - sqrt2a/sqrt1) + (1./w+1./r)*(-sqrt1b/sqrt2 + sqrt2b/sqrt1));
}

double yytFunctions::dgtdd()
{
	return (-2*(sqrt(((1 + d - n)*(1 + m - n/2.))/(1 - m - n/2.)) +
	        sqrt((d*(-m + n/2.))/(m + n/2.)))*
	      (sqrt(((1 + d - n)*(1 - m - n/2.))/(1 + m - n/2.)) +
	        sqrt((d*(m + n/2.))/(-m + n/2.))))/(-2*(-m + n/2.)*(m + n/2.) + n) +
	   ((sqrt(((1 + d - n)*(1 + m - n/2.))/(1 - m - n/2.)) +
	        sqrt((d*(-m + n/2.))/(m + n/2.)))*
	      ((1 - m - n/2.)/
	         (2.*sqrt(((1 + d - n)*(1 - m - n/2.))/(1 + m - n/2.))*(1 + m - n/2.)) +
	        (m + n/2.)/(2.*(-m + n/2.)*sqrt((d*(m + n/2.))/(-m + n/2.))))*(-2*d + n))/
	    (-2*(-m + n/2.)*(m + n/2.) + n) +
	   (((1 + m - n/2.)/
	         (2.*(1 - m - n/2.)*sqrt(((1 + d - n)*(1 + m - n/2.))/(1 - m - n/2.))) +
	        (-m + n/2.)/(2.*sqrt((d*(-m + n/2.))/(m + n/2.))*(m + n/2.)))*
	      (sqrt(((1 + d - n)*(1 - m - n/2.))/(1 + m - n/2.)) +
	        sqrt((d*(m + n/2.))/(-m + n/2.)))*(-2*d + n))/(-2*(-m + n/2.)*(m + n/2.) + n);

//	double sqrt1a = sqrt((1.-w)*(1.-n+d)/(1.-r));
//	double sqrt1b = sqrt(w*d/r);
//	double sqrt1 = sqrt1a + sqrt1b;
//	double sqrt2a = sqrt((1.-r)*(1.-n+d)/(1.-w));
//	double sqrt2b = sqrt(r*d/w);
//	double sqrt2 = sqrt2a + sqrt2b;
//
//	return 0.5*gt()*( dgsdd()/gs() + (sqrt1a/sqrt2+sqrt2a/sqrt1)/(1.-n+d) + (sqrt1b/sqrt2+sqrt2b/sqrt1)/d);
//	printf("++++++++++++ %e, %e\n", aa-c, c/aa);
}


double yytFunctions::dgsdm()
{
	return (-2.*(-2.*(-m + n/2.) + 2.*(m + n/2.))*pow(-2.*d + n,2.))/
			     pow(-2.*(-m + n/2.)*(m + n/2.) + n,3.);
//	return -8*m*gs()/(n-2*r*w);
}

double yytFunctions::dgsdd()
{
	return (-4.*(-2.*d + n))/pow(-2.*(-m + n/2.)*(m + n/2.) + n,2.);
//	return -4*gs()/(n-2*d);
}

double yytFunctions::gamma(double kx, double ky)
{
	return 2.*(cos(kx) + cos(ky));
}

double yytFunctions::eta(double kx, double ky)
{
	return 2.*(cos(kx) - cos(ky));
}

double yytFunctions::epsilon(double kx, double ky)
{
	return -(gt()*t + 3./8.*gs()*J*chi)*gamma(kx,ky);
}

double yytFunctions::BDelta_d()
{
	//printf("in yytF: %f %f %f = %f\n", gs(), J, BDelta, 3/8*gs()*J*BDelta);
	return 3./8.*gs()*J*BDelta;
}

double yytFunctions::BDelta_af()
{
	return 2.*gs()*J*m;
}
